Solving mixed integer programs with Julia/JuMP

Copyright 2016, Pedro Belin Castellucci,

This work is licensed under a Creative Commons Attribution 4.0 International License.</i>

In this Notebook, we will explore some basics of the JuMP modeling language. Although it is a modeling language, it enables the usage of different solvers. By ways of an example, we will use Clp -- a complete list of supported solvers and JuMP documentation can be found at this link. Also, JuMP support linear programming, mixed-integer programming, second-order conic programming, semidefinite programming, and nonlinear programming, but we will focus on linear and integer programming here.

The first step is to install JuMP.


In [ ]:
Pkg.add("JuMP")

We can also install Clp and Cbc in the same way.


In [35]:
Pkg.add("Clp");
Pkg.add("Cbc");


INFO: Nothing to be done
INFO: METADATA is out-of-date — you may not have the latest version of Clp
INFO: Use `Pkg.update()` to get the latest versions of your packages
INFO: No packages to install, update or remove
INFO: Package database updated
INFO: METADATA is out-of-date — you may not have the latest version of Cbc
INFO: Use `Pkg.update()` to get the latest versions of your packages

To get a feeling of JuMP basic functions we will implement the following model:

$Max\ x + y$

subject to:

$2x + 3y \leq 10,$

$3x + 2y \leq 10.$

We begin by creating a Model object.


In [7]:
using JuMP  # Tell that we are using the JuMP library.
using Clp   # Tell that we are using the Clp solver.

m = Model(solver=ClpSolver())  # Creating the Model object.


Out[7]:
$$ \begin{alignat*}{1}\min\quad & 0\\ \text{Subject to} \quad\end{alignat*} $$

The second step is to define the variables. Their are associated with the model and can be created using the <i>@variable</i> macro.


In [8]:
@variable(m, x)
@variable(m, y)


Out[8]:
$$ y $$

Then, the objective function.


In [9]:
@objective(m, Max, x + y)


Out[9]:
$$ x + y $$

Finally, the constraints.


In [10]:
@constraint(m, 2x + 3y <= 10)
@constraint(m, 3x + 2y <= 10)

m  # Just to diplay our model


Out[10]:
$$ \begin{alignat*}{1}\max\quad & x + y\\ \text{Subject to} \quad & 2 x + 3 y \leq 10\\ & 3 x + 2 y \leq 10\\ & x\\ & y\\ \end{alignat*} $$

Now that we have completely defined our model, we may want to solve it.


In [11]:
solve(m)


Out[11]:
:Optimal

The solve function returns the status of the solving procedure. We can query for the solution and variable values using getobjectivevalue and getvalue functions.


In [12]:
println("Objective function = ", getobjectivevalue(m))

println("x = ", getvalue(x))
println("y = ", getvalue(y))


Objective function = 4.0
x = 2.0
y = 2.0

Exercise

Implement and solve the following model:

$Min\ 2x - 3y + z$

subject to:

$x - y - z = 10,$

$3x - 4y + 4z \geq 25,$

$x, y, z \geq 0.$

The least absolute method

In this example, let us assume we have some experimental data and want to fit a linear model to it. Perhaps, the most famous method to accomplish this is the Least Square Method, which finds the linear model that minimizes the squared error, i. e. the distance from data points to the "model points". Instead of minimizing the squared error, we will minimize the absolute error. More formally, we have experimental points $(x_i, y_i)$, $i = 1, \ldots, n$ and we want the parameters $(a, b)$ of a function $f(x) = ax + b$ such that:

$$\sum_{i=1}^n |f(x_i) - y_i|$$

is minimized. We can formulate a linear programming model of the problem using $a$, $b$ and $epsilon$ as decision variables:

$\displaystyle Min\ \sum_{i=1}^n\epsilon_i$,

subject to:

$\displaystyle y_i - a x_i - b \leq \epsilon_i, \quad i=1, \ldots, n$,

$\displaystyle - (y_i - a x_i - b) \leq \epsilon_i, \quad i=1, \ldots, n$ .

To test it, let us generate some data.


In [13]:
n = 10  # We will use 10 points.

points = Array{Float64, Float64}[]
points = zeros(0, 2)

for i=1:n
    
    # They are generated around the curve f(x) = x.
    xRand = rand()*3 - 3
    yRand = rand()*3 - 3

    x, y = i + xRand, i + yRand
    points = [points; x y]
end

points


Out[13]:
10×2 Array{Float64,2}:
 -0.53498  -0.818845
  1.45412   1.88413 
  1.14171   0.368657
  1.13399   1.38762 
  3.01344   4.43634 
  5.28236   4.29358 
  5.43464   5.40538 
  7.82955   6.72188 
  6.24162   8.61881 
  9.169     7.41751 

To get a sense of the data, let us scatter plot it. For that, we are going to use the Plots library.


In [14]:
Pkg.add("Plots");


INFO: Nothing to be done
INFO: METADATA is out-of-date — you may not have the latest version of Plots
INFO: Use `Pkg.update()` to get the latest versions of your packages

And then, we do a scatter plot of the points.


In [15]:
using Plots

plot(points[:, 1], points[:, 2], t=[:scatter], label="", xlabel='x', ylabel='y')


Out[15]:

Now, we implement the model presented earlier. Note how we used the function sum.


In [16]:
using JuMP

m = Model(solver=ClpSolver())  # Creating the model object.

# Defining the variables:
@variable(m, err[1:n])  # error (epsilon).
@variable(m, a)  # angular coefficient.
@variable(m, b)  # linear coefficient.

@objective(m, Min, sum(err[i] for i in 1:n))

for i in 1:n
    @constraint(m, points[i, 2] - a*points[i, 1] - b <= err[i])
    @constraint(m, -(points[i, 2] - a*points[i, 1] - b) <= err[i])
end

m


Out[16]:
$$ \begin{alignat*}{1}\min\quad & err_{1} + err_{2} + err_{3} + err_{4} + err_{5} + err_{6} + err_{7} + err_{8} + err_{9} + err_{10}\\ \text{Subject to} \quad & 0.5349797310141184 a - b - err_{1} \leq 0.8188445915507234\\ & -0.5349797310141184 a + b - err_{1} \leq -0.8188445915507234\\ & -1.454119530901453 a - b - err_{2} \leq -1.8841260173473047\\ & 1.454119530901453 a + b - err_{2} \leq 1.8841260173473047\\ & -1.1417066134268303 a - b - err_{3} \leq -0.36865650956179197\\ & 1.1417066134268303 a + b - err_{3} \leq 0.36865650956179197\\ & -1.1339946573373432 a - b - err_{4} \leq -1.3876237934949982\\ & 1.1339946573373432 a + b - err_{4} \leq 1.3876237934949982\\ & -3.0134434799525147 a - b - err_{5} \leq -4.43633551032433\\ & 3.0134434799525147 a + b - err_{5} \leq 4.43633551032433\\ & -5.282355519486574 a - b - err_{6} \leq -4.293580486595906\\ & 5.282355519486574 a + b - err_{6} \leq 4.293580486595906\\ & -5.434637352223886 a - b - err_{7} \leq -5.405376492945175\\ & 5.434637352223886 a + b - err_{7} \leq 5.405376492945175\\ & -7.829548322209627 a - b - err_{8} \leq -6.721884646395883\\ & 7.829548322209627 a + b - err_{8} \leq 6.721884646395883\\ & -6.2416184000389645 a - b - err_{9} \leq -8.61880516153296\\ & 6.2416184000389645 a + b - err_{9} \leq 8.61880516153296\\ & -9.168997526920457 a - b - err_{10} \leq -7.417506659293617\\ & 9.168997526920457 a + b - err_{10} \leq 7.417506659293617\\ & err_{i} \quad\forall i \in \{1,2,\dots,9,10\}\\ & a\\ & b\\ \end{alignat*} $$

We want to solve the model and get the values of $a$ and $b$.


In [17]:
println("We have obtained an $(solve(m)) solution.")
println("a = $(getvalue(a))")
println("b = $(getvalue(b))")


We have obtained an Optimal solution.
a = 0.7966870433563517
b = 0.4841849427590115

Let us plot the solution to see the result.


In [18]:
fig = plot(points[:, 1], points[:, 2], t=[:scatter], label="")

y = [getvalue(a)*points[i, 1] + getvalue(b) for i in 1:size(points,1)]
plot!(points[:, 1], y, w=4, xlabel="x", ylabel="y", label="")

fig


Out[18]:

Out of curiosity, we can compare our result with a standard Least Square regression. For this, we can use the GLM and DataFrames packages.


In [19]:
Pkg.add("GLM")
Pkg.add("DataFrames")


INFO: Nothing to be done
INFO: METADATA is out-of-date — you may not have the latest version of GLM
INFO: Use `Pkg.update()` to get the latest versions of your packages
INFO: Nothing to be done
INFO: METADATA is out-of-date — you may not have the latest version of DataFrames
INFO: Use `Pkg.update()` to get the latest versions of your packages

In [20]:
using DataFrames
using GLM
data = DataFrame(X=points[:, 1], Y=points[:, 2])
result = glm(Y ~ X, data, Normal(), IdentityLink())
bLS, aLS = coef(result) 

fig = plot(points[:, 1], points[:, 2], t=[:scatter], label="")

y = [getvalue(a)*points[i, 1] + getvalue(b) for i in 1:size(points,1)]
plot!(points[:, 1], y, label="Ours", w=3)

y = [aLS*points[i, 1] + bLS for i in 1:size(points,1)]
plot!(points[:, 1], y, label="LSM", w=3)

fig


WARNING: Method definition describe(AbstractArray) in module StatsBase at /home/pedro/.julia/v0.5/StatsBase/src/scalarstats.jl:560 overwritten in module DataFrames at /home/pedro/.julia/v0.5/DataFrames/src/abstractdataframe/abstractdataframe.jl:407.
WARNING: The bare formula syntax lhs ~ rhs is deprecated. Use @formula(lhs ~ rhs) instead.
 in depwarn(::String, ::Symbol) at ./deprecated.jl:64
 in @~(::Any, ::Any) at /home/pedro/.julia/v0.5/DataFrames/src/statsmodels/formula.jl:32
 in include_string(::String, ::String) at ./loading.jl:441
 in include_string(::Module, ::String, ::String) at /home/pedro/.julia/v0.5/Compat/src/Compat.jl:174
 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/pedro/.julia/v0.5/IJulia/src/execute_request.jl:154
 in #invokelatest#24(::Array{Any,1}, ::Function, ::Function, ::ZMQ.Socket, ::Vararg{Any,N}) at /home/pedro/.julia/v0.5/Compat/src/Compat.jl:505
 in eventloop(::ZMQ.Socket) at /home/pedro/.julia/v0.5/IJulia/src/eventloop.jl:8
 in (::IJulia.##18#24)() at ./task.jl:360
while loading In[20], in expression starting on line 4
Out[20]:

We can see the curves are similar. The least square method have the characteristic of heavily penalize outliers, so points to far from the mean may have a big influence on regression. Note that we have generated the points randomly, you can rerun the experiments to see possible differences in the results.

A more classical example -- The knapsack problem

Let $i \in I$ be an item with value $v_i$ and volume $c_i$. We want to choose the most valuable subset of items to carry in a knaspack without violating its capacity $C$. Let $x_i \in \{0, 1\}$, $i \in I$ indicate whether item $i$ is put into the knapsack. The following integer program solver our knapsack problem.

$\displaystyle Max\ \sum_{i \in I} v_i x_i$

subject to

$\displaystyle \sum_{i \in I} c_i x_i \leq C,$

$x_i \in \{0, 1\}, i \in I.$

We will implement this program using JuMP. First, we need some input data.


In [21]:
v = [3, 1, 4, 6, 2, 7, 2, 4, 7, 9, 5, 7, 3, 1]  # Value of the items.
c = [2, 5, 3, 8, 9, 6, 4, 8, 2, 2, 1, 4, 2, 3]  # Volume of the items.
n = length(v)  # How many items we have

C = 12  # The capacity of our knapsack.


Out[21]:
12

We do the same as before, create our model.


In [37]:
using Cbc
m = Model(solver=CbcSolver())


Out[37]:
$$ \begin{alignat*}{1}\min\quad & 0\\ \text{Subject to} \quad\end{alignat*} $$

As before, we are not going to creating the variables one by one, we will use the following.


In [38]:
@variable(m, x[i=1:n], Bin)


Out[38]:
$$ x_{i} \in \{0,1\} \quad\forall i \in \{1,2,\dots,13,14\} $$

Which is a shortcut for adding the variables x[1] to x[|n|] to the model m. To add the objective function we can use the sum function.


In [39]:
@objective(m, Max, sum(v[i]*x[i] for i=1:n))


Out[39]:
$$ 3 x_{1} + x_{2} + 4 x_{3} + 6 x_{4} + 2 x_{5} + 7 x_{6} + 2 x_{7} + 4 x_{8} + 7 x_{9} + 9 x_{10} + 5 x_{11} + 7 x_{12} + 3 x_{13} + x_{14} $$

It can also be used to create constraints.


In [40]:
@constraint(m, sum(x[i]*c[i] for i=1:n) <= C)

m  # Just to print our model.


Out[40]:
$$ \begin{alignat*}{1}\max\quad & 3 x_{1} + x_{2} + 4 x_{3} + 6 x_{4} + 2 x_{5} + 7 x_{6} + 2 x_{7} + 4 x_{8} + 7 x_{9} + 9 x_{10} + 5 x_{11} + 7 x_{12} + 3 x_{13} + x_{14}\\ \text{Subject to} \quad & 2 x_{1} + 5 x_{2} + 3 x_{3} + 8 x_{4} + 9 x_{5} + 6 x_{6} + 4 x_{7} + 8 x_{8} + 2 x_{9} + 2 x_{10} + x_{11} + 4 x_{12} + 2 x_{13} + 3 x_{14} \leq 12\\ & x_{i} \in \{0,1\} \quad\forall i \in \{1,2,\dots,13,14\}\\ \end{alignat*} $$

And then, we can solve our model.


In [41]:
solve(m)


Out[41]:
:Optimal

Finally, we can query for the solution.


In [42]:
println("Take \t Don't take")

for i in collect(1:n)
    if getvalue(x[i]) >= 0.99
        println("$i")
    else
        println("\t $i")
    end
end


Take 	 Don't take
	 1
	 2
3
	 4
	 5
	 6
	 7
	 8
9
10
11
12
	 13
	 14

Solver callbacks

Some solvers enable a setup of user-defined behavior of some of their internal functionalities. This customization is provided as callback functions. Each solver provides its own set of callback function but JuMP can handle few of them. Here, we will see an example of how to use the Lazy Constraints callback that is supported by Gurobi and CPLEX, among others. For the complete list of callbacks JuMP is able to handle see this link.

As a motivation, we will use the Traveling Salesman Problem (TSP). The TSP consists of finding the least cost route to visit the nodes in a graph exactly once. Let $c_{ij}$ be the distance between nodes $i$ and $j \in V$. Then, binary variables $x_{ij}, i, j \in V, i \neq j$ indicate whether the route uses the edge $(i, j)$. Thus, the objective function of the problem is:

$$Min \ \sum_{i \in V} \sum_{j \in V} c_{ij} x_{ij},$$

in which $c_{ij}$ are the distance (cost) of going from city $i$ to $j$. Two sets of constraints for the problem are:

$$\sum_{i \in V} x_{ij} = 1, \quad j \in V,$$$$\sum_{j \in V} x_{ij} = 1, \quad i \in V.$$

These constraints tell that each node is visited exactly once. We also need to prevent solutions with sub-cycles, this can be done with:

$$\sum_{(i, j) \in C} x_{ij} \leq |C| - 1, \quad C \subset V.$$

Note that we have a constraint for each subset of the nodes, which give a number of constraints that is exponential on the number of nodes. Adding all these constraints a priori to our model could make it too big to be solvable. So we are going to use the Lazy Constraints callback.

To implement it, first, we will generate some data and use the Graphs and Gurobi libraries. To install them, we can use:


In [59]:
Pkg.add("Graphs");
Pkg.add("Gurobi");


INFO: Nothing to be done
INFO: Nothing to be done

Generating the data:


In [58]:
using JuMP
using Gurobi
using Graphs;
# using GraphPlot 

nCities = 9  # Number of cities

dist = zeros(nCities, nCities)

citiesPos = Dict{Int, Any}()

# Position of each city
citiesPos[1] = (0, 0)
citiesPos[2] = (10, 0)
citiesPos[3] = (20, 0)

citiesPos[4] = (0, 10)
citiesPos[5] = (10, 10) 
citiesPos[6] = (20, 10)

citiesPos[7] = (0, 20)
citiesPos[8] = (10, 20)
citiesPos[9] = (20, 20)

for (i, posI) in citiesPos
    for (j, posJ) in citiesPos
        dist[i, j] = ((posI[1] - posJ[1])^2 + (posI[2] - posJ[2])^2)^0.5
    end 
end

# Plotting the nodes of the graph:
g = simple_graph(nCities, is_directed=true)
# Graphs.plot(g)


could not spawn `neato -Tx11`: no such file or directory (ENOENT)

 in _jl_spawn(::String, ::Array{String,1}, ::Ptr{Void}, ::Base.Process, ::Base.PipeEndpoint, ::Base.DevNullStream, ::Ptr{Void}) at ./process.jl:321
 in #424 at ./process.jl:478 [inlined]
 in setup_stdio(::Base.##424#425{Cmd,Ptr{Void},Base.Process}, ::Tuple{Pipe,Base.DevNullStream,Base.PipeEndpoint}) at ./process.jl:466
 in #spawn#423(::Nullable{Base.ProcessChain}, ::Function, ::Cmd, ::Tuple{Pipe,Base.DevNullStream,Base.PipeEndpoint}, ::Bool, ::Bool) at ./process.jl:477
 in (::Base.#kw##spawn)(::Array{Any,1}, ::Base.#spawn, ::Cmd, ::Tuple{Pipe,Base.DevNullStream,Base.PipeEndpoint}, ::Bool, ::Bool) at ./<missing>:0
 in open(::Cmd, ::String, ::Base.DevNullStream) at ./process.jl:544
 in #plot#43(::String, ::Function, ::Graphs.GenericGraph{Int64,Graphs.Edge{Int64},UnitRange{Int64},Array{Graphs.Edge{Int64},1},Array{Array{Graphs.Edge{Int64},1},1}}) at /home/pedro/.julia/v0.5/Graphs/src/dot.jl:99
 in plot(::Graphs.GenericGraph{Int64,Graphs.Edge{Int64},UnitRange{Int64},Array{Graphs.Edge{Int64},1},Array{Array{Graphs.Edge{Int64},1},1}}) at /home/pedro/.julia/v0.5/Graphs/src/dot.jl:91
 in include_string(::String, ::String) at ./loading.jl:441

Now, we implement the method that will tell whether a solution found by the solver has a sub-cycle. The function must receive a parameter which is the callback handler. Then, we get the solution found and try to find its sub-cycles. We also must use the <i>@lazyconstraint</i> macro to actually add the constraint to the model and re-optimize.


In [45]:
function subtourElimConstraint(cb)

    # Creating a graph:
    g = simple_graph(nCities, is_directed=true)

    for i in 1:nCities, j in 1:nCities
        if i != j
            if getvalue(x[i, j]) > 0.99
                add_edge!(g, i, j)  # Populating the graph.
            end
        end
    end
    
    # Plot it, so we can see the graph:
    Graphs.plot(g)
    
    nodes = vertices(g)
    
    # The idea is to walk through the edges and check how many edges
    # we have to hop to reach the origin. This number will tell us if
    # we are walking on a sub-cycle. If that is the case, we add the 
    # corresponding constraint to the model.
    # We repeat this procedure starting at an unvisited node for finding
    # multiple sub-cycles.
    
    visited = []
    for n in nodes

        if n in visited
            continue
        end

        pivot = n  # Starting node
        out = out_edges(pivot, g)[1]
        outNode = target(out, g)
        append!(visited, [pivot, outNode])

        subtourConstraint = AffExpr()
        subtourConstraint += x[pivot, outNode]
        tourSize = 1

        # While we don't get to the origin, keep hoping...
        while(pivot != outNode)
            out = out_edges(outNode, g)[1]
            newOut = target(out, g)

            subtourConstraint += x[outNode, newOut]
            tourSize += 1
            push!(visited, outNode)
            outNode = newOut
        end

        # If we have traversed the whole graph:
        if tourSize == num_vertices(g)
            return
        else
            # This is the constraint we will add:
            println(subtourConstraint, " <= ", tourSize - 1)
        end
        
        # Using the @lazyconstraint macro
        @lazyconstraint(cb, subtourConstraint <= tourSize - 1)        
    end
    
end


Out[45]:
subtourElimConstraint (generic function with 1 method)

We have the method that will add constraints as we need them. Let us check how this works on the model implementation. We begin as usual, defining model, variables, the objective function and constraints...


In [ ]:
model = Model(solver=GurobiSolver(TimeLimit=20, Threads=1))

@variable(model, x[i=1:nCities, j=1:nCities; i != j], Bin)

@objective(model, Min, sum(dist[i, j] * x[i, j] for i in 1:nCities, j in 1:nCities if i != j));

for i in 1:nCities
    @constraint(model, sum(x[i, j] for j in 1:nCities if i != j) == 1)
end 
        
for j in 1:nCities
    @constraint(model, sum(x[i, j] for i in 1:nCities if i != j) == 1)
end

Then, before solving, we need to tell that we will be using the callback and which is the function with the customized behavior.


In [ ]:
addlazycallback(model, subtourElimConstraint)

# And the, solving:
solve(model)

To check the final solution, we can proceed as follows.


In [ ]:
# Creating a graph:
g = simple_graph(nCities, is_directed=true)

for i in 1:nCities, j in 1:nCities
    if i != j
        if getvalue(x[i, j]) > 0.99
            add_edge!(g, i, j)  # Populating the graph.
        end
    end
end

# Plot it, so we can see the graph:
Graphs.plot(g)

Note that the graph library we are using do not plot the nodes in their exact positions. You can draw them to check if it makes sense.

In this Notebook, we have explored the basics of JuMP regarding mixed integer linear programming and the lazy constraint callback. For a more detailed exposition the reader is referred to JuMP documentation.

Hope you had fun! (=